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We examine avalanche statistics of rain- and vibration- driven granular slides 
in miniature sand mounds. A crossover from power-law to non power-law avalanche- 
size statistics is demonstrated as a generic driving rate v is increased. For 
slowly-driven mounds, the tail of the avalanche-size distribution is a power- 
law with exponent —1.97 ± 0.31, reasonably close to the value previously 
reported for landslide volumes. The interevent occurrence times are also an- 
alyzed for slowly-driven mounds; its distribution exhibits a power-law with 
exponent -2.670 ± 0.001. 
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13 Introduction. — Landslides are the movement of a mass of rock, debris, or earth down a 

u slope triggered by a variety of natural factors, ranging from rainfall to volcanic activity 

15 On 17 February 2006, a series of mudslides caused widespread damage and loss of life in 

16 Southern Leyte, Philippines. The deadly landslides followed a ten-day period of persistent 

17 downpour and a minor (M2.6) earthquake [Catane et ai, 2007]. 

18 One practical approach in analyzing the underlying physical processes that generate 

19 landslide statistics is numerical modeling. One class of models, based on self- organized 

20 criticality (SOC) [Hergarten and Neugebauer, 1998; Piegari et ai, 2006], hint at some 

21 possible mechanisms yielding the observed statistics. SOC is a theory underlying the 

22 spontaneous emergence of critical-like behavior (i.e., power laws and critical exponents) 

23 in systems for which the timescales between buildup and release of stress are separated, 

24 and for which the stress-transfer mechanism is generally nonconservative [Juanico et ai, 

25 2007a, b; Juanico and Monterola, 2007]. SOC concepts have aroused great interest in 

26 the study of granular matter [Jaeger et al, 1989], a well-known example of which is the 

27 ricepile experiment [Frette et ai, 1996]. 

28 The present work examines avalanche statistics of rain- and vibration- driven granular 

29 slides in miniature sand mounds. A previous study [Katz and Aharonov, 2006] explored 

30 slope-failure types due to horizontal and vertical vibrations in a miniature sandbox. It 

31 was shown that vertical shaking leads to a power-law distribution of slide-block surface 

32 area, although the experiment did not demonstrate a 'rollover' observed in substantially 

33 complete, empirical landslide inventories [Guzzetti et al., 2002; Malamud et al., 2004a, b]. 

34 In the present study, we incorporate rainfall as a triggering mechanism and show that 
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its presence allows our model to capture the rollover. Likewise, we demonstrate ex- 
perimentally the existence of a crossover from power-law to non power-law statistics as 
theoretically predicted by Piegari et al. [2006]. 

It appears from the data reported by Malamud et al. [2004a] that any type of trigger 
(earthquake, rainfall, or snowmelt) yields roughly the same trend in the landslide-size 
distribution. Combining any two trigger types should reasonably yield the same trend in 
the distribution, and this is precisely one of the aspects tested in this study. In addition, 
the computational model proposed here assumes a generic triggering mechanism. Thus, 
rainfall and vibrations have been introduced to act as concurrent landslide triggers. 

Experimental Method. — The experimental setup (Figure 1) consists of a sand mound 
(total mass, 1500 g) disturbed concurrently by rainfall and vibration. River sand with 
irregularly- shaped grains (mean grain mass, 8.8 x 10~ 4 g; mean grain volume, 5.8 x 
10~ 4 cm 3 ) was used. Initially, the mound is dry and is a near-perfect cone in shape with a 
base diameter of 23 cm and height of 8 cm (slope angle ~ 35° with respect to horizontal). 
Water (200 ml) is dispensed quite uniformly over the surface and allowed time to seep 
into the mound. The wet mound thus consists of 12% water by weight at the start of the 
observations. 

Rainfall is simulated using a makeshift sprinkler (water capacity, 75 ml) placed directly 
above the mound apex. Water pours out through a circular orifice of diameters: 4.5, 10.0, 
and 47.0 (±0.5) mm. Pour rate is the mass of water flowing out per unit time, and the 
orifice diameter controls its value. Pour rate is constant at the following values: 3.21 g s _1 ; 
6.60 g s _1 ; and 21.09 g s" 1 over a time interval of 20, 10, and 3 s, respectively. These 
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57 time intervals are commensurate with the length of time the 75-ml water in the sprinkler 

58 IS depleted. 

59 Horizontal vibrations are applied using a tabletop earthquake simulator. Horizontal 
eo shaking at low to moderate accelerations generates grain flows of the topmost layer of the 

61 slope resulting in rapid failure-plane development [Katz and Aharonov, 2006]. As shown 

62 in Figure 1, a translational-load platform, powered by a servo motor driven by a USB 6009 

63 DAQ driver (National Instruments™), imparts the horizontal vibrations having a saw- 

64 tooth wave profile. The wave profile is fed into the DAQ driver via LabVIEW™computer 
es interface. The wave is characterized by a maximum amplitude of 1.5 cm (with respect 
66 to center) and by frequencies of: 1.8 Hz, 10 Hz, and 89 Hz. Due to the 3D shape of the 
e? mound, slope-parallel and slope-normal accelerations are both present during shaking. 

68 Computational Model. — The underlying physics of the avalanche statistics of driven 

69 granular slides is investigated by performing numerical experiments of a landslide model 

70 proposed by Piegari et al. [2006], defined as follows. A mountain 'slope' is represented as 

71 a 2D inclined plane partitioned into a grid (500 cells x 500 cells). Each cell k is described 

72 by a stress parameter 9k initialized randomly between and 1 from an arbitrarily chosen 

73 rectangular distribution (although the distribution used for initial randomization does 

74 not affect the long-term behavior of the model [Piegari et a/., 2006]). The randomization 

75 captures the expected heterogeneity of stress values in actual mountain slopes. Stress in 

76 the slope builds up over time by means of a localized (i.e., cell scale) driving: 9k (t + At) = 

77 9f.(t) + v At, where v is the generic driving rate. For simplicity, it is assumed that v has 

78 the same value for all cells. When a cell k has 9^ > 1, it relaxes by transferring stress to 
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its four nearest neighbors nn = {up, down, right, left} at different proportions g nn . By 
virtue of gravity, stress transfer is biased downwards, such that: gd own > 9up, subject to 
the constraint gdown + 9u P = 0.5; and gi e f t = g r i g ht = 0.25. In this study, it is assumed that 
transfer is conservative, so that ^Z nn g n n = 1- Stress-transfer proceeds until the entire 
grid relaxes so that for all k, 9^ < 1. All consecutive stress-transfers comprise a landslide 
at time t, and the total number of collapsing cells at time t is the landslide area A(t). 
For correspondence with experiments, although we recognize its limitations, we assume a 
landslide mass-area relation: M ~ A 3 / 2 , which can be derived by means of dimensional 
analysis (i.e., assuming that M is proportional to volume V, and then considering that 
V ~ A 3/2 [Hovius et al, 1997]). 

The novelty of our computational modeling approach is in considering the gradual 
flattening of the slope after several landslides have occurred (as seen in our experiments). 
An update rule is introduced to decrease gdown, as follows: gdownit + At) = gdownit) — 
10 _5 A(t), where A(t) is the area of landslide at time t. The value of g up is updated 
accordingly via the constraint gdown + g up = 0.5. This modification incorporates the 
dynamics of slope evolution which was not realized by Piegari et al. [2006]. The generic 
driving rate v defined in the model corresponds to the experimental parameters, as shown 
in Table 1. 

Results and Discussion. — Avalanche size is interpreted as the total mass M of wet sand 
falling from the mound onto the basin within a 20-ms period (temporal resolution of our 
data-capturing device). M is measured by a weighing scale (resolution, 0.1 g) beneath 
the basin, as illustrated in Figure 1. M is sometimes contributed to by several distinct 



DRAFT 



August 25, 2008, 12:19am 



DRAFT 



JUANICO ET AL.: AVALANCHE STATISTICS OF DRIVEN GRANULAR SLIDES... X - 7 

101 'sub-avalanches.' The tradeoff in this interpretation is that although it resolves mostly 

102 individual avalanches, those that last > 20 ms may be recorded as partial sub-avalanches, 

103 and if there are more than one avalanche within any 20 -ms interval, these will be recorded 

104 as one value. 

105 The total observation period is determined by the pour rate and is based on how long 

106 before the entire mound washes out. A total of 20, 10, and 3 trials were made for slow, 
10? moderate, and fast pour rates, respectively. The fact that sand tends to stick together 
loa when wet (i.e., negative pore pressure) [Katz and Aharonov, 2006] justifies our assumption 
109 that by measuring the mass M of falling wet sand, the avalanche volume is effectively mea- 
uo sured (i.e., M oc V). In [Malamud et al, 2004a, b], the distribution of landslide volumes 
in has been deduced from scaling arguments due to the difficulty in obtaining direct informa- 
n2 tion about landslide volume from field data [Malamud et al, 2004a]. In our experiment, 
us landslide volume distribution is determined by using mass as a proxy for volume. 

n4 Figure 2 illustrates probability densities (pdf) of avalanche sizes resulting from slow (♦), 

us moderate (•), and fast (■) driving (as defined in Table 1). Also shown are corresponding 

He pdfs from numerical simulations of the computational model. The pdfs for all cases exhibit 

in peaks that shift towards the right as the driving rate is increased. In particular, The tail of 

us the pdf for a slowly-driven mound (♦, Fig. 2) is a power-law with exponent — 1.97(±0.31; 

U9 standard error of least-squares fit, df = 4, reduced \ 2 — 1-14 x 10~ 5 ), determined by 

120 curve-fitting the linear portion (in double log scale) for which 1.0 g < M < 8.0 g. The 

121 exponent value is not significantly different (two-tailed i-test: t = 1.13, df = 4, p = 0.32) 
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122 from the value —1.93 reported for probability densities of landslide volumes [Malamud 
et ai, 2004a]. 

124 A crossover from power-law to non power-law statistics is also evident in our data. This 

125 crossover has been previously recognized by Piegari et al. [2006] as an effect of increasing 

126 the generic driving rate v. We attribute the crossover to timescale separation. Based 
12? on the model, the timescale separation is the ratio between the timescale v~ x of stress 
us changes on any site due to driving and the timescale of the longest avalanche which is set 

129 at At = 1. Hence, the timescale separation is (z/At)" 1 ~ v~ x . For slow driving (small v\ 

130 SOC theory expects the appearance of power laws. On the other hand, for fast driving, 
m a different trend is expected. We found this fit to resemble a normal distribution, and is 

132 centered at a large M (■, Fig. 2). An implication of the crossover is that under high-rate 

133 driving (especially by frequent rainfall), the landslide behavior of mountain slopes appears 

134 to produce relatively large avalanches on average. It is thus not surprising that highly- 

135 devastating rain-induced landslides occur more often in areas frequently struck by storms. 
o6 However, a more comprehensive description of v should incorporate ground-failure factors 
is? such as soil composition. Such factors have been neglected in this study for simplicity. 

us The rollover observed in our data supports the claim that the rollover seen in sub- 

139 stantially complete landslide-inventory datasets is not a mere artifact of limited mapping 

«o resolution [Guzzetti et ai, 2002; Malamud et ai, 2004a, b]. As v — >• , the rollover is 

i4i expected to disappear thus leaving the power-law tail of the distribution. This agrees 

«2 with predictions from the Olami-Feder-Christensen model [Olami et ai, 1992], which is 

143 the limiting case of our model for v — > 0. To the extent that our experiment and numerical 
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«4 results agree, we thus attribute the rollover to the finiteness, albeit small, value of the 

145 driving rate v assumed to characterize real physical systems. From this argument it can 

we thus be reasonably deduced that the rollover increases in prominence, such that the peak 

14? shifts towards larger avalanche sizes, as v increases. We found that at v ~ 7.5 x 10 -4 , 

us the avalanche size distribution more closely resembles a normal (or Gaussian) distribution 

149 that is centered at a large size value. 

150 A relevant aspect for hazard prediction is the interevent occurrence time (IOT) statis- 

151 tics. The IOT is the interval between the peaks of events whose sizes are above a given 

152 threshold (M = 0.5 g), which corresponds roughly to the peak of the pdf (♦, Figure 2). 

153 We thus gathered time series of avalanche size M(t) in slowly-driven mounds over a period 

154 of 20 s. Figure 3 illustrates a representative sample for the first 10 s of this time series for 

155 both experiment and computational data. A thresholding procedure is applied to discard 
ise events whose sizes are below the threshold. The region M < 0.5 g corresponds to the 
is? rollover portion of the pdf (♦, Figure 2) for slowly-driven mounds. The vertical demar- 
iss cation line in Figure 2 corresponds to the horizontal demarcation line in Figure 3. The 
159 probability density for interevent occurrence time derived from experiment and model are 
wo shown on the inset graph of Figure 3 (A; blue, solid curve), and both exhibit a power-law 

161 trend with exponent — 2.670(±0.001; standard error of least-squares fit, df = 6, reduced 

162 x 2 = 2.09 x 10~ 8 ). 

163 A power-law IOT distribution implies that most correlated events (i.e., those with sizes 

164 above the threshold) tend to occur close together in time, which seems to agree with 

165 recent findings [Rossi et al, 2008]. The power-law trend further implies that correlated 
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lee events may be separated a long time from each other — a landslide today may be linked 

16? with an earlier landslide a long time ago. This may be tied to the fact that landslides tend 

lea to occur where they have occurred before. The long temporal correlation suggested by 

169 the power- law could possibly be attributed to rainfall seasonal trends [Rossi et al, 2008]. 

no To illustrate that a power-law emerges from temporal correlations in the time series 

171 data, a shuffling procedure has been implemented wherein the order of the time series 

172 is rearranged randomly [Yang et al, 2004]. Random shuffling effectively eliminates any 

173 trace of correlations present in the original time series data. The thresholding procedure 

174 is then applied on the shuffled data to extract the IOT distribution. After shuffling, the 

175 IOT distribution becomes exponential, as shown in Figure 3 (▼, Inset). The change from 

176 power-law to exponential due to shuffling has been expected for systems governed by SOC 

177 dynamics [Woodward et al, 2004]. Therefore, we conclude that the power-law is a direct 

178 indication of temporal correlations in the time series data. 

179 Conclusion. — The general agreement of our results with empirical data suggests that 
lac miniature experimental models may help in understanding several underlying facets of 

181 complex landslide processes. While several geophysical factors such as sand porosity, 

182 rock type, pore-water pressure, and humidity have been neglected, the landslide model 

183 presented here delivers basic insights that would guide more detailed explorations later 
™ on. 
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Table 1. Generic driving rate v in terms of experimental parameters 



Description 


Driving rate (no units) 


V 


Pour rate 




Vibration frequency (Hz) 


Slow 


7.50 x 10' 


-5 




3.21 


1.8 


Moderate 


1.25 x 10" 


-4 




6.60 


1.8, 10, 89 


Fast 


7.50 x 10" 


-4 




21.09 


1.8, 10, 89 



DRAFT 



August 25, 2008, 12:19am 



DRAFT 



X - 14 



JUANICO ET AL.: AVALANCHE STATISTICS OF DRIVEN GRANULAR SLIDES. 



«1 



□ □□□ □ n □ □ □ □ □ 
Desktop □□□□ □□□□ nnn 
nn ™ a DI.Il 



Ml 







\ 1 


Basin 


1 1 °-°°a 1 r> 


Weighing 
Scale 


Servo Motor 



Nl USB 6009 DAQ Power amplifier 



Translational 
Load Stage 



Figure 1. Diagram of experimental setup. Sand mound consists of river sand. Rainfall 
is simulated by a sprinkler right above the mound apex. An earthquake is simulated by a 
translational load stage imparting horizontal vibrations to the platform where the mound 
is placed. Rainfall and vibration are parametrized by pour rate and vibration frequency, 
respectively. Different parameter combinations are listed on Table 1. Avalanche size is 
the mass M of wet sand falling onto the basin within 20-ms intervals. M is measured 
by a weighing scale at the bottom of the basin, and is connected to the PC that records 
measurements in time. 
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Figure 2. (Color online). Probability densities (pdf) of avalanche sizes. Pdfs for slow 
(♦), moderate (•), and fast (■) driving; where the error bars are ±2a/t2/(5 (equivalent to 
±2cr, 95% confidence; n = number of values in a bin; and § is the bin width). Corre- 
sponding pdfs of avalanche size from numerical simulations (curves) are overlaid, after 
converting area into mass using the scaling relation: M = (1.0 x 10 _5 g ■ cel\~ 3 ^ 2 )A 3 ^ 2 . 
Linear regime (1.0 < M < 8.0 g) of the pdf for a slowly-driven mound fits a power-law 
with exponent — 1.97 ± 0.31. For moderately-driven mounds, the rollover region of the 
pdf becomes more prominent as the peak shifts to the right. For highly-driven mounds, 
the pdf resembles a normal distribution centered at a large M. 
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Time (s) 



Figure 3. (Color online). Interevent occurrence times (IOT). Time series of avalanche 
sizes from experiment (red, solid curve) and model (blue, thick curve). Consecutive data 
points are separated by 20 ms. Iteration steps (iter) in the model have been rescaled to 
time by multiplying by the scaling factor: 0.05 s- iter -1 . Avalanche sizes less than 0.5 g are 
discarded. IOT is the interval between consecutive peaks above the threshold. Inset: IOT 
distributions for experiment (a) and model (blue solid curve) both exhibit a power-law 
with exponent —2.670 ± 0.001. Randomly shuffling the order of the time series effectively 
results in an exponential IOT distribution (▼; curve fit: green, dashed curve). 
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